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Following the recent realization of an artificial version of Graphene in the electronic surface states of copper 
with judiciously placed carbon monoxide molecules inducing the honeycomb lattice symmetry (K. K. Gomes 
et a/.Nature 483, 306 (2012)), we demonstrate that these can be used to realize a vortex in a Kekule texture of 
the honeycomb lattice. The Kekule texture is mathematically analogous to a superconducting order parameter, 
opening a spectral gap in the massless Dirac point spectrum of the Graphene structure. The core of a vortex 
in the texture order parameter, supports subgap states, which for this system are analogs of Majorana fermions 
in some superconducting states. In particular, the electron charge bound to a single vortex core is effectively 
fractionalized to a charge of e/2. The Kekule texture as realized in the molecular Graphene system realizes 3 
different domain types, and we show that a Y-junction between them realizes the coveted Kekule vortex. 

PACS numbers: 



The experimental realization of Graphene 1 , has inspired in 
the last few years a large body of work exploring the possi- 
ble physics in an ideal 2D system, realizing massless Dirac 
fermions, including a great deal of theoretical work 2 . While 
Graphene is proving a very flexible medium to manipulate, as 
a physical system it has its limitations, and sadly, some of the 
most interesting physical effects theoretical work suggested 
in Graphene-like systems have not been realized in the actual 
Graphene system. For this reason, there is good reason to ex- 
plore alternative realizations of the single layer honeycomb 
2D electron gas. 

Recent advances in STM technology have allowed to man- 
ufacture an artificial version Graphene, by arranging car- 
bon monoxide molecules on the surface of copper, dubbed 
"molecular Graphene". The carbon monoxide molecules are 
arranged in a regular array, and thus create an electrostatic 
potential with minima forming the vertices of a honeycomb 
lattice*^ 4 -. Any 2D electron gas with the symmetry of the 
honeycomb lattice imposed on it is likely to realize an ana- 
log of Graphene. Because of its microscopic construction, 
the molecular Graphene system is even more easily manip- 
ulated than Graphene. In particular, since the electrostatic 
potential is essentially under full control by selecting an ap- 
propriate molecule arrangement, the honeycomb lattice can 
be engineered with a wide variety of lattice textures, which 
are predicted to realize the analog of huge magnetic fields^—, 
as well as analogs of superconducting states* --^. Among the 
most intriguing of these proposals is the suggestion to realize 
a Kekule texture on the honeycomb lattic e) 10 - 14 : 15 . The Kekule 
texture makes some nearest neighbor links on the honeycomb 
lattice stronger than others, in a x s/3 arrangement as de- 
picted in Fig. [TJ In the low energy effective massless Dirac 
Hamiltonian for the honeycomb tight binding model, this ar- 
rangement induces a pairing-like term between the electrons 
of one Dirac point and the holes of the second Dirac point (in- 
stead of electrons with opposite spin). In principle the Kekule 
order parameter can have a complex value, and can include 
a vortex. It has been shown* , that this vortex is the analog 
of a vortex in a p x + ip y superconducting state, supporting a 



zero mode at the vortex core***. In analogy to some systems 
proposed in the context of high energy physic a 17 - 18 , and other 
solid state system o 19 ' 20 the vortex core is expected to bind a 
fermion number of one half, in the Kekule texture case - half 
an electron, if spin is ignored. The halving of a fermion is 
also at the heart of the emergence of Majorana fermions in 
superconducting states - there the fermion being halved is the 
Bogoliubov quasiparticle. The halved electron realized by the 
Kekule vortex would be a mathematical analog of a Majorana 
fermion. In this paper we demonstrate how a Kekule vortex 
can be realized in the molecular Graphene system, opening 
up the possibility to explore directly the physics of fermion 
halving. 

At the microscopic level a uniform Kekule texture enlarges 
the unit cell of the honeycomb lattice 3 fold (as illustrated in 
Fig- ID- The unit cell includes three plaquettes, one of which 
has the nearest neighbor hopping strength on the links around 
it stronger (or weaker). With any choice of unit cell there are 
three choices where to locate the Kekule distorted plaquettes. 
We will show that the Y-junction between these three domains 
realizes a vortex in the Kekule texture, binding a charge e/2 
(per spin) to its core. 

Even without the texture, we can use the same enlarged unit 
cell to describe the Graphene band structure, from which we 
can identify the low energy massless Dirac Hamiltonian. Af- 
terwards, we add the Kekule texture. 

The enlarged unit cell, as depicted in Fig. [2] has 6 lattice 
sites, instead of 2 in the original unit cell. We therefore de- 
note the six sites in the unit cell by /i = 1 ... 6, at positions 
a p = (sin (^p) , cos (^fp)) relative to the unit cell center. 
The original Bravais lattice basis of ai + a 2 and a,5 + 8l$ is 
replaced with Ai = 3ai and A2 = 3a2. Finally, the en- 
larged real space unit cell, corresponds in momentum space 
to copying the band structure with a shift of the reciprocal 

lattice vectors B x = ^ i^' 1 ) and B 2 = x (t!' -1 )' 
The first Brillouin Zone is 3 times smaller, and both Dirac 
points are shifted to q = 0, since the Dirac points appear at 
Q = ±(Bi +B 2 ). 
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FIG. 1: Kekule pattern on the honeycomb lattice. The thick (red) 
links have stronger (or weaker) hopping strength on them. The unit 
cell is denoted by dashed (gray) lines. 




FIG. 2: The unit cell of the Kekule pattern. The six sites in the 
unit cell are denoted by 1 ... 6 (red). The unit cell is denoted by the 
dashed (gray) line. The (blue) Bravais lattice vectors are denoted by 
Ai 2 . 



The tight binding Hamiltonian is 



n ° = t J2J2 c U r ) |>n( r ) + s-i( r ) + s+3( r + A ^ 

r 



(i) 



where r denotes the unit cell, A M = 3a M are Bravais lat- 
tice vectors, and the /i ± 1, fj, + 3 indices are added mod- 
ulo 6. The hopping strength is set to t = 1 for conve- 
nience. Going to momentum space via the Fourier transform 
c„(r) = J e? 2 qe 4( ' r+a ^- ) ' q c M (q) , the Hamiltonian becomes 

Ho = J d\ cl(q)H^(q)cM) > (2) 



where 



tfrf(q) 



(3) 



The eigenstates are found from the equation iJ (q)^(q) = 

e(q)V>(q)- 

The key to identifying the Dirac points in this matrix 
form is to perform an appropriate basis change. We know 
that we could have chosen a smaller unit cell with just 2 
sites (wavefunction weights Xi,2), and a choice of ^(q) = 

(Xi, X2, Xie 1 ** 3 - 1 , X2e 1 ^ , Xie 1 ** 4 - 1 , X 2 e 4qa6 - 2 ) T , 
where q = 0, ±Bi, and ay = a; — &j should recover 
that choice, since it recovers the plane wave phases be- 
tween the different sites in the unit cell, while keeping the 
2 sites of the small unit cell unchanged. This therefore 
suggests using the unitary transformation comprised of 
U = (Vi(0) ) V 2 (0) ) ^i(Bi) j ^(B 1 ),Vi(-Bi) j V2(-Bi)) 
, where Vi(q) = (l, 0, e^* 3 "* 1 ), 0, e^ C* 4 "* 1 ), 0) T , and 
-0 2 (q) = (0, l,0,e iq '( a4 - a2 \0,e lc i-( a6 - a2 )). This turns out 
to be 
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where to = e i27r / 3 . Using this unitary transformation, expand- 
ing to linear order in q and taking t = I, we find 
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where q± = q± ± zg2- The structure of the two Dirac points 
is now easily seen in the diagonal blocks in the third to sixth 
columns. The diagonal 2x2 block in the first and second 
columns corresponds to high energy modes we will ignore for 
the low energy theory. 

Next we add the Kekule texture to the Hamiltonian, with 
strength A. With any choice of unit cell there are three choices 
of the Kekule pattern, shown in Fig. [3] With our choice of 
unit cell the additional hopping strength in each case can be 
quantified by adding to Ho the terms if,\,a=i,2,3, where 
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iq-(a„-a M ) (S^^S^oda + ^/i-l.i/^.even) + e !q ' a "^+3^ 

(6) 
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corresponding to the patterns in Fig. |3(a)| Fig. |3(b)| and 
Fig. |3(c)| respectively. 

Taking now H = Hq + H\ for the three cases, expand- 
ing to linear order in small q, and finally applying the unitary 
transformation (0), we find for H\i 



£/t (H + H Xil )U = 






2A 





|wg+ 





-|w 2 g 


2 A 










-|w 2 q+ 








±uj 2 q+ 





M- 





—Aw 







Xq + 





—Aw 













-Aw 2 





-Ag+ 







-Aw 2 





-A<?_ 






(7) 

where A = 1 + A. We project out the high energy modes 
involving the first and second columns and rows, and retain 
only the low energy Hamiltonian blocks, which yield 

(H + H x , a =i,2, 3 )U = 
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We find that the Kekule texture renormalizes the Fermi veloc- 
ity by A = 1 + A, and otherwise gives a term mixing between 
the two Dirac cones. The mixing term has a different phase 
for the three domains, giving a Kekule texture order parame- 
ter A = Aw a , where a = 1, 2, 3. The relative differences in 
phase being From this fact we conclude that at a Y junc- 
tion between the three different domains of Kekule texture, 
there will be a phase winding of 2ir, thus realizing a vortex. 

Next we present a numerical calculation of the local density 
of states (LDOS) near the Kekule vortex. The LDOS is a par- 
ticularly useful quantity in the context of molecular Graphene 
since STM is used to construct the molecular Graphene to 
begin with, and the same STM can be used to measure the 
LDOS. We will use this to try and probe the unique states 
bound to the Kekule vortex. We use the tight binding model 
Ho on a finite, disc-shaped flake of Graphene (of radius 22.5 
in the convention we use here), with Kekule textures added 
realizing the Y-junction between the three domains of the 
Kekule texture, as shown in Fig. [4] Note that a disc geom- 
etry was chosen to minimize spurious states appearing at the 
system edges (for instance at corners). We calculate the eigen- 
states of the system W\ip a ) — e a \^Pa), with A = 1, and find 
the LDOS using the formula 

v»(E,r) =^5(E - e a )\(0\ Cli (T)\1i a )\ 2 , (9) 

a 

where S(x) is in the ideal case is a Dirac delta function, but 
for a calculation in a finite system must be taken as some 
approximation for the delta function. We take a Lorentzian 

$( x ) — x 2 "+w 2 °f width w = 0.00001 as our approximation. 

We plot the LDOS on the various lattice sites for a number 
of different energies E in Fig. [5] and find that for E — the 
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FIG. 3: The three different Kekule domains, with a specific choice 
of unit cell. These three domains are characterized by an effective 
phase difference of 2n/3. 




FIG. 4: Image of the Kekule texture Y-junction, which we use in 
numerics. The Graphene flake is a disc of radius 22.5, using the 
convention mentioned in the text. The stronger hopping links are 
denoted by a thick (red) link between the (gray) points denoting the 
lattice sites. The boundaries between the three Kekule texture do- 
mains are denoted by (black) lines, indicating the Y-junction shape. 
The point where the three domain walls meet, is a vortex core. An 
effective phase jumps by 2tv/3 at each domain wall. 



LDOS is strongly peaked at the vortex center, and in a spot at 
the edge of the system (see Fig. |5(d)| each realizing (roughly) 
one half of an electron. This is precisely what one expects in 
the case of a halved fermion in a finite system. In the ideal 
case, a zero mode appears bound to the vortex core, and an- 
other zero mode appears bound to the system edge. These 
states are degenerate in energy, and any infinitesimal matrix 
element between them will cause them to form symmetric 
and anti-symmetric linear combinations, slightly split in en- 
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ergy. The lower of these two energy states is a fermion state 
delocalized between two positions, with half its wavefunction 
weight at each spot, regardless of the distance between them. 

Next we will want to find numerically the charge accumu- 
lated at the vortex core. We do this in two ways. Taking the 
lower energy E « state,we can integrate its wavefunction 
weight up to some radius R away from the vortex core 

D 1 (R)= [ R rdr [ +X d^ -(r)\ 2 = E l^-(rj)| 2 , 

J ° \r s \<R 

(10) 

where r, <fr are the polar coordinates in the plane. Also, we can 
take the LDOS at E = 0, and integrate it up to R 

D 2 (R) = / rdr dcf>v{T,E = 0) = ^ v{r j ,E = Q). 
Jo J - 7T \ rj \<R 

(11) 

It is important to note that the latter is a quantity we have 
experimental access to. We plot the calculated D\ t % (R) versus 
R in Fig. [6] Ideally, only the zero modes at the vortex core 
and the edge will contribute to the LDOS at E = 0. In a finite 
system, we will get the sum of the contributions from the two 
slightly split linear combinations of the zero modes 

i/(r,£7 = 0)« / dEv{v,E) = |^ _(r)| 2 + |^ 0+ (r)| 2 . 
Jo- 

(12) 

It is therefore expected that half the weight of this quantity 
will be at the vortex core, and the other half at the edge. While 
Di(R) saturates at a value of 1, D2{R) saturates at an arbi- 
trary value, and indicates that we will not get a quantitative 
measure of the accumulated charge in the vortex core from 
the LDOS. However, we can still learn a great deal from 
the quantity STM can measure in the lab, as it does demon- 
strate that roughly half the overall weight is accumulated in 
the vortex core, the remaining weight being concentrated near 
the disc edge. Observing roughly half the total weight cen- 
tered at the vortex core would suggest that a fractionalized 
state exists, but this is not conclusive. 

Experimental measurement is further complicated by the 
fact that the electron spin needs to be taken into account. 
In the carbon monoxide on copper system for molecular 
Graphene, spin orbit coupling is negligible, and interactions 
seem to be weak, and therefore all electronic states ought to 
be spin degenerate, we will still have electron halving be- 
tween a vortex core and the sample edge, but this will occur 
for both spin polarizations. The LDOS measured by STM 
would be the sum of the contributions from the two spins, but 
we should still observe a curve like that of Fig. |6(b)| Further- 
more, we can move the spin up and down states in opposite 
direction by applying a Zeeman field, sufficiently weak not to 



cross any other electronic state, but sufficiently strong to split 
the different spin zero mode states sufficiently to be observed 
in the LDOS measurement. 

An additional complication arises from the fact that the real 
system also has a non vanishing second neighbor hopping t' 
which lowers the symmetry of the Hamiltonian (|T), and ruins 
the theoretically perfect e/2 fractionalization 10,21 , changing it 
to some more general fraction. However, the qualitative dis- 
tinction of the vortex core bound states remains - one electron 
is delocalized between the vortex core and the edge of the sys- 
tem, with some finite fraction of its weight bound to the vortex 
core, and the the rest to the edge. 

Perhaps a better method than merely measuring the static 
LDOS, averaged over long times, would be to probe some 
correlation between the edge and the vortex core, or bet- 
ter yet between two vortex cores. As explained earlier, the 
effective fermion halving is essentially the delocalizing of 
a single fermion wavefunction between two spots (for in- 
stance two vortex cores). Qualitatively, if the electron is de- 
tected near one of the vortex cores at some short time inter- 
val, then the wavefunction collapses onto that vortex core, 
and no electron should be detected at the other vortex core 
during the same short time interval. An experiment prob- 
ing this temporal correlation could perhaps reveal the funda- 
mental quantum mechanical effect at play here. One could 
try to simultaneously measure time resolved electron tunnel- 
ing at two locations. Calculating the noise correlation be- 
tween the two tunneling currents 1 1,2, averaged over time 
(Aii A/2) = (Iih) — (Ii)(l2), should reveal a long range 
correlation, only between the vortex cores. Detecting this 
would be a strong indicator that a wavefunction at E ps 
is delocalized between these two locations, thus realizing the 
fermion halving scenario. 

In conclusion we have demonstrated that the molecular 
Graphene system can be made to form a Kekule texture with 
a vortex, thus realizing a physical system with fermion halv- 
ing. In this case the electron effectively fractionalizes to states 
with charge e/2 bound to the vortex core. The electron spin is 
expected to merely double the electronic spectrum, and thus 
the vortex core should accumulate a unit charge, but no mag- 
netization due to spin (see also Ref.lToh. The Kekule texture 
Y-junction has already been realized experimentally^ 2 ., and it 
now remains to prove that a fermion halving is indeed occur- 
ring in this system. We explored how signatures of the halv- 
ing would appear in the LDOS, and hope our insights will be 
tested in the molecular Graphene system. 
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FIG. 5: The LDOS for A = 1, at different energies: e = +1.5[(a)| 
e = +0.5[(b)] e = +0.1^ e = 0[(d)] e = -0.l[(e)] e = -0.5 
[(f)] and e = — 1.5 |(g)| The circles represent sites of the honeycomb 
lattice, taken here in the shape of a disc of radius 22.5, using the con- 
ventions in the main text. The site coloring is such that dark (blue) 
points have a higher weight, and lighter (orange) points have lower 
weight. While LDOS scans above [(a)! and below [(g)! the gap show a 
rather uniform distribution of DOS, in the gap there clearly is some 
spatial structure. In particular, at e = we find a peak at the vortex 
core, and at one spot on the disc edge |(d)[ as expected. 
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FIG. 6: Radial accumulated weight of the zero mode wavefunction 
Di(fl)[(a)l and of the LDOS at E = D 2 (R ) [(b)] The probability 
density of the V'o- wavefunction is depicted in |(c)| using the same 
convention as for the LDOS plots in Fig. [5] 



